#include <stdio.h>
#ifdef USEMPI
#include<mpi.h>
#endif
#ifdef USEOMP
#include<omp.h>
#endif
long int  num_steps = 1000000000;
double step;
int main (int argc, char* argv[])
{   
int id=0,numprocs=1,provided;
long int i;
double x, pi, sum = 0.0;
double r_sum =0.0;   //for mpi only
long int low, up , res;

step = 1.0/(double) num_steps;
res  = num_steps%numprocs;

#ifdef USEMPI
    MPI_Status status;      
    MPI_Init_thread(&argc,&argv,MPI_THREAD_FUNNELED, &provided); 
    MPI_Comm_rank(MPI_COMM_WORLD,&id);
    MPI_Comm_size(MPI_COMM_WORLD,&numprocs);
#endif 
if (res == 0) {
    low = id*(num_steps/numprocs)+1 ;
    up  =  low + num_steps/numprocs-1;
   } else 
   {
    if(id == 0) { 
        low = 1; 
        up = num_steps/numprocs;}
    else if(id <= res) {      
       low = id*(num_steps/numprocs)+id; 
       up = low + num_steps/numprocs; }
    else if(id > res) {
       low = id*(num_steps/numprocs)+1+res; 
       up = low + num_steps/numprocs-1;
       } 
   }

//printf(" 1 pi is %2.20f\n",pi);  
#ifdef USEOMP 
#pragma omp parallel for reduction(+:sum) private(i,x)
#endif
    for (i=low;i<= up; i++){   
        x = ((double)i-0.5)*step;
        sum = sum + 4.0/(1.0+x*x);
        }
    sum=sum*step;
    r_sum = sum;
#ifdef USEMPI
    MPI_Barrier(MPI_COMM_WORLD);
    MPI_Reduce(&sum,&r_sum,1, MPI_DOUBLE, MPI_SUM,0,MPI_COMM_WORLD) ; 
    MPI_Barrier(MPI_COMM_WORLD);
#endif
    if(id == 0) {  
        pi = r_sum;
       printf("pi is %2.20f\n",pi);
    }
#ifdef USEMPI  
    MPI_Finalize();
#endif 
    return 0;
}
